# Librarieslibrary(tidyverse)library(gmRi)library(rnaturalearth)library(scales)library(sf)library(gt)library(patchwork)library(lmerTest)library(emmeans)library(merTools)library(tidyquant)library(ggeffects)library(performance)library(gtsummary)library(gt)library(sizeSpectra)library(ggdist)library(pander)library(ggh4x)conflicted::conflict_prefer("select", "dplyr")conflicted::conflict_prefer("filter", "dplyr")conflicted::conflicts_prefer(lmerTest::lmer)# Processing functionssource(here::here("Code/R/Processing_Functions.R"))# ggplot themetheme_set(theme_gmri(rect =element_rect(fill ="white", color =NA), strip.text.y =element_text(angle =0),axis.text.x =element_text(size =8),legend.position ="bottom"))# vectors for factor levelsarea_levels <-c("GoM", "GB", "SNE", "MAB")area_levels_long <-c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")area_levels_all <-c("Northeast Shelf", area_levels)area_levels_long_all <-c("Northeast Shelf", area_levels_long)# table to join for swapping shorthand for long-hand namesarea_df <-data.frame(area =c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),survey_area =c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),area_titles =c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))# Degree symboldeg_c <-"\u00b0C"# Function to get the Predictions# Drop effect fits that are non-significant ###get_preds_and_trendsignif <-function(mod_x){ modx_preds <-as.data.frame(# Model predictionsggpredict( mod_x, terms =~ est_year * survey_area * season) ) %>%mutate(survey_area =factor(group, levels = area_levels),season =factor(facet, levels =c("Spring", "Fall")))# Just survey area and year modx_emtrend <-emtrends(object = mod_x,specs =~ survey_area*season,var ="est_year") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Preds with signif modx_preds %>%left_join(select(modx_emtrend, survey_area, season, non_zero))}
Code
# Broad Distribution#log2 bins for easy-of-access#' @title Build Log 2 Bin Structure Dataframe#' #' @description Used to build a dataframe containing equally spaced log2 bins for#' size spectra analysis. Contains details on the left and right limits, midpoint, bin width, #' and a text label for the bins. log2bin number ascends with increasing size for eeasy plotting.#'#' @param log2_min #' @param log2_max #' @param log2_increment #'#' @return#' @export#'#' @examplesdefine_log2_bins <-function(log2_min =0, log2_max =13, log2_increment =1){# How many bins n_bins <-length(seq(log2_max, log2_min + log2_increment, by =-log2_increment))# Build Equally spaced log2 bin df log2_bin_structure <-data.frame("log2_bins"=as.character(seq(n_bins, 1, by =-1)),"left_lim"=seq(log2_max - log2_increment, log2_min, by =-log2_increment),"right_lim"=seq(log2_max, log2_min + log2_increment, by =-log2_increment)) %>%mutate(bin_label =str_c(round(2^left_lim, 3), " - ", round(2^right_lim, 3), "g"),bin_width =2^right_lim -2^left_lim,bin_midpoint = (2^right_lim +2^left_lim) /2) %>%arrange(left_lim)return(log2_bin_structure)}#' @title Assign Manual log2 Bodymass Bins - By weight#'#' @description Manually assign log2 bins based on individual length-weight bodymass #' in increments of 1 on the log2 scale. Returns data with bins assigned based on individual#' length-weight biomass#' #' Uses maximum weight, and its corresponding bin as the limit.#'#' @param wmin_grams Catch data prepared for mle calculation, use prep_wmin_wmax#' @param log2_increment Equally spaced increments to use for log 2 bin sizes. Default = 0.5.#'#' @return#' @export#'#' @examplesassign_log2_bins <-function(wmin_grams, log2_increment =1){#### 1. Set up bodymass bins# filter missing weights size_data <- wmin_grams %>%filter(wmin_g >0,is.na(wmin_g) ==FALSE, wmax_g >0,is.na(wmax_g) ==FALSE)# Get bodymass on log2() scale size_data$log2_weight <-log2(size_data$ind_weight_g)# Set up the bins - Pull min and max weights from data available#min_bin <- floor(min(size_data$log2_weight)) min_bin <-0 max_bin <-ceiling(max(size_data$log2_weight))# Build a bin key, could be used to clean up the incremental assignment or for apply style functions log2_bin_structure <-define_log2_bins(log2_min = min_bin, log2_max = max_bin, log2_increment = log2_increment)# Loop through bins to assign the bin details to the data log2_assigned <- log2_bin_structure %>%split(.$log2_bins) %>%map_dfr(function(log2_bin){# limits and labels l_lim <- log2_bin$left_lim r_lim <- log2_bin$right_lim bin_num <-as.character(log2_bin$log2_bin)# assign the label to the appropriate bodymasses size_data %>%mutate(log2_bins =ifelse( between(log2_weight, l_lim, r_lim), bin_num, NA),log2_bins =as.character(log2_bins)) %>%drop_na(log2_bins) })# Join in the size bins log2_assigned <-left_join(log2_assigned, log2_bin_structure, by ="log2_bins")# return the data with the bins assignedreturn(log2_assigned)}
Storycrafting - Northeast Shelf Spectra
I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills
This document will be split into two sections (the second potentially moving to its own doc).
Part 1: Summary of Trends
For the first section the focus is on understanding changes in the size structure. This will be done using length and bodymass spectra, large and small fish indices, and looking at how body mass size quantiles have moved over time.
This section should address: Has the size distribution changed over time, where has it changed, and some nuance around what aspects about the distribution have shifted.
Starting from Square Uno
I’d like for us to start by rebuilding the story of how the size structure has changed over time and regional differences.
As such, I would like for us to look at plots of the size spectrum slope for (1) the all-species length spectrum and (2) the Wigley species biomass spectrum from annual, seasonal and regional perspectives.
I am thinking this would be two panel plots (all-spp length, Wigley-spp biomass) with 5 subpanels each (NES, GOM, GB, SNE, MAB) that are set up to show lines by seasonal and annual time steps, much like the plot on the right of slide 7 in your most recent slide deck…except with NES added and an annual line added. Keeping the significant trend lines in is helpful. This is the first step, so when you have that ready, please share it with the group. - KMills
Code
# Data used for Wigley estimateswigley_in <-read_csv(here::here("Data/processed/wigley_species_trawl_data.csv"))finfish_in <-read_csv(here::here("Data/processed/finfish_trawl_data.csv"))
Aside: Biomass Fraction in Each Community
If we need to make a decision about which community we want to use, here is the different biomass totals for each:
Code
# glimpse(wigley_in)# Load raw# General tidying only, removal of strata outside our study area, Spring and Fall only# (removes inshore and rarely/inconsistently sampled strata etc.)# shrimps removedtrawl_basic <-read_csv(here::here("Data/processed/trawl_clean_all_data.csv"))# Biomass totalsraw_totals <- trawl_basic %>%distinct(est_towdate, cruise6, station, stratum, tow, comname, biomass_kg) %>%pull(biomass_kg) %>%sum()wig_biomass_total <- wigley_in %>%distinct(est_towdate, cruise6, station, stratum, tow, comname, biomass_kg) %>%pull(biomass_kg) %>%sum()ffish_biomass_total <- finfish_in %>%distinct(est_towdate, cruise6, station, stratum, tow, comname, biomass_kg) %>%pull(biomass_kg) %>%sum()# Make a table for displaytibble("Community"=c("All (including crustaceans)", "All finfish", "Wigley Species Subset"),"Biomass Total"=c(raw_totals, ffish_biomass_total, wig_biomass_total)) %>%mutate(`% of Total`= (`Biomass Total`/ raw_totals)*100) %>%gt()
Community
Biomass Total
% of Total
All (including crustaceans)
3781548
100.00000
All finfish
3708046
98.05629
Wigley Species Subset
3539719
93.60501
A. Size Spectra Trends
For the figures in part 1: non-zero trends over time have been highlighted in the appropriate plot panels with linear trends overlaid.
# Left side:# All species species length spectrum # color = season vs. annual# facet_wrap(~survey_area)# Right side:# Wigley species biomass spectrum# color = season vs. annual# facet_wrap(~survey_area)# Plot over observed data# Contrast seasonal differences( (lenspectra_trend_p +theme(strip.text.y =element_blank()) | bmspectra_trend_p) ) +plot_layout(guides ="collect")
B. Median Length & Weight Trends
I also think it would be valuable to build out the story of size change further, so am thinking that plots like those of average (or median) length (all species) and weight (wigley species) for NES and the sub-regions, like those on p. 12 of the attached file above (the very first draft from late 2022) would be useful. - KMills
Code
# Load the median weight/length datawigley_size_df <-read_csv(here::here("Data/processed/wigley_species_size_summary.csv"))finfish_size_df <-read_csv(here::here("Data/processed/finfish_species_length_summary.csv"))# shelf-wide summaryfinfish_lengths_shelf <-read_csv(here::here("Data/processed/shelfwide_finfish_species_length_summary.csv"))wigley_sizes_shelf <-read_csv(here::here("Data/processed/shelfwide_wigley_species_size_summary.csv"))# Join themsize_results <-bind_rows(list("Finfish Community"=bind_rows(finfish_size_df, finfish_lengths_shelf),"Wigley Subset"=bind_rows(wigley_size_df, wigley_sizes_shelf)), .id ="community") %>%mutate(survey_area =factor(survey_area, levels = area_levels_all) )# Get trends:len_mod_ffish <-lmer(formula = med_len_cm ~ est_year * survey_area * season + (1|est_year),data = finfish_size_df)# len_mod_wig <- lmer(# formula = med_len_cm ~ est_year * survey_area * season + (1|est_year),# data = wigley_size_df)wt_mod_wig <-lmer(formula = med_wt_kg ~ est_year * survey_area * season + (1|est_year),data =drop_na(wigley_size_df, med_wt_kg))# Get the predictions and flag whether they are significantsize_trend_predictions <-bind_rows(list("med_len_cm-Finfish Community"=get_preds_and_trendsignif(len_mod_ffish),"med_wt_kg-Wigley Subset"=get_preds_and_trendsignif(wt_mod_wig) ), .id ="model_id") %>%separate(model_id, into =c("metric", "community"), sep ="-") %>%mutate(metric =fct_rev(metric))# Left side: # Median length - all species# color = season vs. annual# facet_wrap(~survey_area)size_long <- size_results %>%pivot_longer(cols =c(med_len_cm, med_wt_kg), names_to ="metric", values_to ="val") %>%mutate(survey_area =fct_relevel(survey_area, area_levels_all),season =fct_rev(season))# just length for scaleslen_plot <- size_long %>%filter(metric =="med_len_cm", community =="Finfish Community") %>%drop_na(val) %>%ggplot() +geom_ribbon(data =filter(size_trend_predictions, metric =="med_len_cm", non_zero ==TRUE),aes(x, ymin = conf.low, ymax = conf.high, fill = season),linewidth =1, alpha =0.35) +geom_point(#data = filter(size_long, metric == "med_len_cm"),aes(est_year, val, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, val, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +geom_line(data =filter(size_trend_predictions, metric =="med_len_cm", non_zero ==TRUE),aes(x, predicted, color = season, linetype ="Significant Trend"),linewidth =1, alpha =1) +scale_fill_gmri() +scale_color_gmri() +facet_grid(survey_area ~ metric*community,scales ="free") +theme(strip.text.y =element_blank()) +labs(x ="Year",y ="Length (cm)",color ="Season",linetype ="Linetype",fill ="Season")# Right: # Median weight - wigley species# color = season vs. annual# facet_wrap(~survey_area)# weight plotswt_plot <- size_long %>%filter(metric =="med_wt_kg") %>%drop_na(val) %>%ggplot() +geom_ribbon(data =filter(size_trend_predictions, metric =="med_wt_kg", non_zero ==TRUE),aes(x, ymin = conf.low, ymax = conf.high, fill = season),linewidth =1, alpha =0.35) +geom_point(aes(est_year, val, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, val, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +geom_line(data =filter(size_trend_predictions, metric =="med_wt_kg", non_zero ==TRUE),aes(x, predicted, color = season, linetype ="Significant Trend"),linewidth =1, alpha =1) +scale_fill_gmri() +scale_color_gmri() +facet_grid(survey_area ~ metric*community,scales ="free") +labs(x ="Year",y ="Weight (kg)",color ="Season",linetype ="Linetype",fill ="Season")(len_plot | wt_plot) +plot_layout(guides ="collect")
There are a number of years when median weight surges in MAB. For these situations I think we likely could trace these down to exceptionally high catches of larger species like elasmobranchs.
We also don’t see a shelf-wide decline in size so I might need to check what Friedland did or prepare for a confrontation.
C. Small/Large Fish Trends (True Percentiles)
Data for part C. comes from the Wigly species subset to set percentiles based on individual bodymass.
The following figures frame “large fish” as individuals with bodymass greater than the 90th percentile for all years within that region. “Small fish” are considered smaller than the 10th percentile of the mass distribution.
The following table shows what those percentiles are for each region.
NOTE: Percentiles are determined using DescTools::Quantile() which uses numlen_adj to weight everything by abundance.
Code
# # Get 90th or 95th percentile as size threshold# Do them separatelyarea_size_quants <-bind_rows(mutate(wigley_in, survey_area ="Northeast Shelf"), wigley_in) %>%split(.$survey_area) %>%map_dfr(function(x,y){ DescTools::Quantile(# x$length_cm, x$ind_weight_kg,weights = x$numlen_adj, probs =c(0.05, 0.1, 0.5, 0.90, .95)) %>%t() %>%as_tibble() }, .id ="survey_area") %>%mutate(survey_area =factor(survey_area, levels = area_levels_all)) %>%arrange(survey_area)# Show the tablearea_size_quants %>%mutate_if(is.numeric, round, 2) %>%gt() %>%tab_header(title ="Wigley Species - Body Mass Quantiles (kg)") %>%tab_spanner(label ="Percentile", columns =-survey_area) %>%cols_label("survey_area"="")
There’s really just so few “large” individuals that the percentiles sit quite small.
Code
# Weighted histograms are possible, we can overlay the median size as wellshelf_median <- wigley_in %>%#filter(ind_weight_g >1) %>% mutate(survey_area ="Northeast Shelf") %>%group_size_metrics(size_data = ., .abund_col ="numlen_adj", .group_cols =c("survey_area", "season"),has_weights =TRUE, .weight_col ="ind_weight_kg")# Plot the weighted histogramwigley_in %>%filter(ind_weight_kg <=10) %>%ggplot(aes(x = ind_weight_g, weight = numlen_adj, fill = season)) +geom_histogram(fill ="#535353") +facet_wrap(~fct_rev(season), nrow =2) +scale_y_continuous(transform =transform_log(base =2),labels =label_log(base =2),expand =expansion(add =c(0,0))) +scale_x_continuous(transform =transform_log(base =2),labels =label_log(base =2),breaks =2^seq(0,16, 1),limits =c(1, 2^15),expand =expansion(add =c(0,0))) +geom_vline(data = shelf_median,aes(xintercept = med_wt_kg*1000, color = season),linewidth =1) +labs(title ="Overall Bodymass Distribution for Wigley Species",subtitle ="Weights over 10kg removed, count not* normalized (expect slope ~0)",x ="Mass (g)",color ="Median")
Code
# We need to normalize to be comparable to the Edwards method# Might be able to just assign the bins and their widths, then divide the weights by that# ^ seems to work!# add the spectra slope values from MLE method as a side value, could use facet grid label# plot the normalizedwigley_in %>%mutate(survey_area =factor(survey_area, levels = area_levels)) %>%filter(ind_weight_kg <=30, est_year %in%c(1990:2019)) %>%assign_log2_bins(log2_increment =1) %>%mutate(season = season) %>%ggplot(aes(x = ind_weight_g, weight = numlen_adj/bin_width)) +geom_histogram(aes(fill =after_stat(density)),breaks =2^seq(0,16,1), color ="gray90", linewidth =0.2) +# geom_vline(# data = filter(wigley_size_df, est_year %in% c(1990:2019),# survey_area != "Northeast Shelf"),# aes(xintercept = med_wt_kg*1000,# color = "Median Weight")) +scale_y_continuous(transform =transform_log(base =2),labels =label_log(base =2),expand =expansion(add =c(0,0))) +scale_x_continuous(transform =transform_log(base =2),labels =label_log(base =2),breaks =2^seq(0,16, 2),limits =c(1, 2^14),expand =expansion(add =c(0,0))) +scale_fill_distiller(palette ="YlGnBu", direction =1) +scale_color_gmri() +facet_nested(est_year~survey_area*fct_rev(season)) +#facet_grid(est_year~survey_area) +theme(axis.text.y =element_blank(), panel.grid.major.y =element_blank(), legend.position ="right") +labs(title ="Normalized Biomass Spectra",x ="Mass (g)",y ="Density",color ="")
The axes on these are: \[\frac{\sum{individuals_{ijs}
> mass~threshold}}
{\sum{individuals_{ijs}}} * 100\] for each: \(year_i\), \(region_j\), \(season_s\).
Code
# Plot them as loing-term average# Plot the SFIsfi_p <- lfi %>%ggplot() +geom_point(aes(est_year, SFI, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, SFI, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_percent()) +theme(strip.text.y =element_blank()) +facet_grid(survey_area~., scales ="free") +labs(x ="Year",y ="Percent of Total Abundance",title ="Small Fish Index\n(< 10th Percentile Bodymass)",fill ="Season",color ="Season",linetype ="")# Plot the LFIlfi_p <-ggplot(lfi) +geom_point(aes(est_year, LFI, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, LFI, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_percent()) +facet_grid(survey_area~., scales ="free") +labs(x ="Year",y ="",title ="Large Fish Index\n(> 90th Percentile Bodymass)",fill ="Season",color ="Season",linetype ="")# Plot together(sfi_p | lfi_p) +plot_layout(guides ="collect")
This plots the abundance of individuals above/below the size thresholds as abundance.
Code
# Plot them as loing-term average# Plot the SFIsfi_p <- lfi %>%ggplot() +geom_point(aes(est_year, small_abund, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, small_abund, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_number()) +theme(strip.text.y =element_blank()) +facet_grid(survey_area~., scales ="free") +labs(x ="Year",y ="Total Abundance",title ="Small Fish Index\n(< 10th Percentile Bodymass)",fill ="Season",color ="Season",linetype ="")# Plot the LFIlfi_p <-ggplot(lfi) +geom_point(aes(est_year, large_abund, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, large_abund, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_number()) +facet_grid(survey_area~., scales ="free") +labs(x ="Year",y ="",title ="Large Fish Index\n(> 90th Percentile Bodymass)",fill ="Season",color ="Season",linetype ="")(sfi_p | lfi_p) +plot_layout(guides ="collect")
D. 1kg Discrete Size Breaks
Due to the skewed distribution in abundance with increasing size, using percentile cutoffs leads to a situation where “large” individuals are in practice quite small.
For the plots below I took the the 1-5kg range in body weights and broke that range into five equally spaced body size ranges.
Code
# What percent is 5kg and smaller?wigley_smol <-filter(wigley_in, ind_weight_kg <=5)smol_frac <-signif((sum(wigley_smol$numlen_adj, na.rm = T) /sum(wigley_in$numlen_adj, na.rm = T)), 3) *100
This still includes smol_frac% of the abundance for wigley species.
Code
# Do everything together# Make evenly spaced binsoverall_weight_ranges <-seq(floor(min(wigley_smol$ind_weight_kg, na.rm = T)),ceiling(max(wigley_smol$ind_weight_kg, na.rm = T)),length.out =6)# Make that a table for plottingoverall_range_tbl <-tibble("size_group"=c("S", "S-M", "M", "M-L", "L"),"min"= overall_weight_ranges[1:5],"max"= overall_weight_ranges[2:6] ) %>%mutate(size_range =str_c(min, "-",max, "kg"))# Plot the breaksoverall_range_tbl %>%mutate(size_group =factor(size_group, levels =c("S", "S-M", "M", "M-L", "L"))) %>%ggplot() +geom_segment(aes(x = min, xend = max, y =0, yend =0, color = size_group), linewidth =14) +geom_point(aes(x = min, y =0), shape ="|", size =5) +geom_point(aes(x = max, y =0), shape ="|", size =5) +scale_y_continuous(expand =expansion(add =c(0.001, 0.001)), labels =NULL, breaks =NULL) + rcartocolor::scale_color_carto_d(palette ="OrYel") +labs(x ="Individual Body Mass (kg)", y ="", color ="Size Range", title ="Even Breaks in Body Mass",subtitle =">5kg Individuals Removed")
If we plot these 1kg bins we can see that 0-1kg individuals dominate the region.
These figures really should be standardized into some sort of CPUE to account for differences in effort and area.
Code
# Plot the Northeast Shelf changeseven_breaks_summary %>%mutate(size_range =fct_rev(size_range)) %>%filter(survey_area !="Northeast Shelf") %>%mutate(area_title ="Northeast Shelf") %>%ggplot(aes(est_year, total_size_abund, fill = survey_area)) +geom_col(position ="stack") +scale_fill_gmri() +scale_y_continuous(breaks = scales::breaks_pretty(), labels =label_number(big.mark =",") ) +facet_grid(size_range*size_group~area_title*season, scales ="free") +labs(y ="Abundance", x ="Year", fill ="Region")
You can really see the jump at 2008 in these figures if you look at Spring, and at the smaller bins for Fall.
Biomass totals below are from biomass-at-length estimates*
Code
# Plot the Northeast Shelf changeseven_breaks_summary %>%mutate(size_range =fct_rev(size_range)) %>%filter(survey_area !="Northeast Shelf") %>%mutate(area_title ="Northeast Shelf") %>%ggplot(aes(est_year, total_size_bio, fill = survey_area)) +geom_col(position ="stack") +scale_fill_gmri() +scale_y_continuous(breaks = scales::breaks_pretty(), labels =label_number(big.mark =",") ) +facet_grid(size_range*size_group~area_title*season, scales ="free") +labs(y ="Biomass (kg)", x ="Year", fill ="Region",title ="Biomass Totals, not effort adjusted")
Summary of Trends
Let’s start with these pieces, and we can then discuss doing something with the extremes of the size distribution–changes in small vs. large individuals, where we pre-define size bins as you’ve done in slides 7 and 8, or perhaps using an approach similar to Lora’s approach of defining the upper/lower percentiles of size. - KMillz
Part 2: Relationships to External Factors
The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.
I will lead with a characterization of broad regional changes in temperature and fisheries landings totals.
Code
# # Model Dataset for Wigley Species Subsetwigley_bmspectra_model_df <-read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))# load shelfwide stuff toowigley_bmspectra_shelf <-read_csv(here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))bot_temps <-read_csv(here::here("Data/processed", "trawl_region_seasonal_bottom_temps.csv"))# Bolt that onwigley_bmspectra_model_df <-bind_rows( wigley_bmspectra_model_df,left_join(wigley_bmspectra_shelf, bot_temps, by =join_by("est_year"=="year", survey_area, season)))# Use Wigley bmass spectra model data for plotswigley_bmspectra_model_df <- wigley_bmspectra_model_df %>%mutate(survey_area =case_when( survey_area =="GoM"~"Gulf of Maine", survey_area =="GB"~"Georges Bank", survey_area =="SNE"~"Southern New England", survey_area =="MAB"~"Mid-Atlantic Bight", survey_area =="Northeast Shelf"~"Northeast Shelf"),survey_area =factor(survey_area, levels = area_levels_long_all),season =fct_rev(season))
Bottom Temperature
Bottom Temperature Trends
Bottom temperature data is from the DuPontavice and Saba combined bottom temperature product. This data product has been clipped to the study areas and seasonally averaged to match the survey sampling seasons.
Plotted below are the long-term seasonal changes. Significant seasonal trends are shown with linear trends and error displayed.
Seasonal bottom temperatures have increased overall for the Northeast Shelf. In all four regions ans shelf-wide there is a significant trend in annual averages, but seasonally spring temperatures in SNE and MAB are stable.
Seasonal Bottom Temperature Distribution
The Gulf of Maine experiences a much smaller difference in seasonal temperature swings than the other regions. The difference in bottom temperature from Spring to Fall can be smaller than the variability of spring bottom temperature along the whole Northeast Shelf.
Code
# Plot the distribution as a raincloud plotwigley_bmspectra_model_df %>%ggplot(aes(fct_rev(survey_area), bot_temp, color = season, fill = season)) + ggdist::stat_halfeye(alpha =0.4,adjust = .5, width = .6, .width =0, justification =-.2, point_colour =NA ) +geom_boxplot(fill ="transparent",width = .15, outlier.shape =NA ) +## add dot plots from {ggdist} package ggdist::stat_dots(## orientation to the leftside ="left", size =1,## move geom to the leftjustification =1.12, ## adjust grouping (binning) of observations binwidth = .1) +scale_color_gmri() +scale_fill_gmri() +scale_y_continuous(labels =label_number(suffix ="\u00b0C")) +coord_flip() +theme(legend.position ="none") +labs(y ="Bottom Temperature", x ="", color ="Season",title ="Seasonal Bottom Temperature Distributions")
The overlap in bottom temperature distributions is very similar to the overlap in mass spectra exponents. Gulf of Maine spectra and bottom temperatures share a good region of overlap, and the regions that experience larger seasonal fluctuations in BT see larger seasonal fluctuations in their size distributions.
Code
wigley_b_dist_plot <- wigley_bmspectra_model_df %>%ggplot(aes(fct_rev(survey_area), b, color = season, fill = season)) + ggdist::stat_halfeye(alpha =0.4,adjust = .5, width = .6, .width =0, justification =-.2, point_colour =NA) +geom_boxplot(fill ="transparent",width = .15, outlier.shape =NA) +## add dot plots from {ggdist} package ggdist::stat_dots(## orientation to the leftside ="left", size =1,## move geom to the leftjustification =1.12, ## adjust grouping (binning) of observations binwidth = .005) +scale_color_gmri() +scale_fill_gmri() +coord_flip() +labs(y ="Exponent of Size Distribution", x ="",title ="Body Mass Distributions",subtitle ="Bodymass Ranges (min: 1g, max: max observed)")# Plot next to bottom temperatureswigley_b_dist_plot
Seasonal Catchability Differences
Bottom temperature relationships to the individual size distribution are emergent at the seasonal time-scale, and broad North-South geographic scales. The extent which the size distribution follows temperature is largely a reflection in the seasonal change in the size spectra. Within any area and for the shelf as a whole, the relationship is much weaker within one season.
This makes me lean towards belief that the size distribution is more responsive on shorter time scales, adjusting on the fly to the ambient water temperatures in order to align so well with seasonal temperature swings. This would also suggest that as a whole the community size spectra may be less sensistive to the long-term trends since the distribution is able to effectively track larger absolute swings in temperature that happen within the year.
I think there is an important question about temperature’s role, but I don’t think that local warming is having much impact on the community size distribution. If it is, it is secondary to other size distribution shifts happening at shorter time scales and across spatial scales.
This then raises the questions of how seasonal size distribution changes are being accomplished. From the decadal variability work we have reason to believe that individual species are both not tracking temperature changes in space at the rate they are ocurring AND that changes in size at age are happening, signals at the species level that temperature effects on growth are ocurring. There would need to be some compensatory mechanisms present to offset this across the broader community.
Total Landings
The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.
The following model has been used to determine the relationship between landings and seasonal bottom temperature on b.
The Northeast Shelf as a whole is not included as its own region in these models, but is plotted below with points just for visual context.
lmer(
b ~ survey_area*season*scale(roll_temp) + survey_area * scale(log(total_weight_lb)) + (1|est_year),
data = wtb_model_df)
Code
#### Model 2: Temperature & Fishing ##### # Drop NA's# wtb_model_df <- # drop_na(wigley_bmspectra_model_df, total_weight_lb, bot_temp, b) %>%# rename(area = survey_area) %>% # left_join(area_df) %>% wtb_model_df <-drop_na(wigley_bmspectra_model_df, bot_temp, b) %>%# Do rolling BT within a season * regiongroup_by(survey_area, season) %>%mutate(roll_temp = zoo::rollapply(bot_temp, 5, mean, na.rm = T, align ="right", fill =NA)) %>%ungroup() %>%mutate(yr_num =as.numeric(est_year),yr_fac =factor(est_year),survey_area =factor(survey_area, levels = area_levels_long_all),season =factor(season, levels =c("Spring", "Fall")),landings = total_weight_lb,yr_seas =str_c(season, est_year))# Make the modelmod2 <-lmer( b ~ survey_area*season*scale(roll_temp) + survey_area*scale(log(total_weight_lb)) + (1|est_year), data = wtb_model_df)# summary(mod2)# check_model(mod2)
Temperature Relationship
Looking at all the data there is a strong relationship between seasonal size distribution and temperature.
Code
btemp_overall <-ggplot(wtb_model_df, aes(bot_temp, b)) +geom_point(aes(shape = survey_area, color = season)) +scale_color_gmri() +geom_smooth(method ="lm", color ="black") +labs(shape ="Area", color ="Season",title ="Overall Temperature Relationship")btemp_overall
But
If we use yearly region * season combinations as our units of comparison the relationship between seasonal bottom temperature and the size distribution falls apart.
Code
# Clean up the plot:# Plot the predictions over datamod2_preds <-as.data.frame(ggpredict( mod2, terms =~ roll_temp*survey_area*season)) %>%mutate(survey_area =factor(group, levels = area_levels_long_all),season =factor(facet, levels =c("Spring", "Fall")))#### Trend Posthoc ####trend_df <-emtrends(object = mod2,~survey_area * season,var ="roll_temp")# Drop effect fits that are non-significant ###mod2_emtrend <-emtrends(object = mod2,specs =~ survey_area*season,var ="roll_temp") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Limit temperature plotting rangeactual_range <- wtb_model_df %>%group_by(season, survey_area) %>%summarise(min_temp =min(bot_temp)-2,max_temp =max(bot_temp)+2,.groups ="drop")# Plot over observed databtemp_local <- mod2_preds %>%left_join(actual_range) %>%filter((x < min_temp) == F, (x > max_temp) == F) %>%left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%filter(non_zero) %>%ggplot() +geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = survey_area), alpha =0.1) +geom_line(aes(x, predicted, color = season), linewidth =1) +geom_point(data = wtb_model_df,aes(roll_temp, b, color = season),alpha =0.4,size =1) +facet_grid(survey_area~., scales ="free") +scale_color_gmri() +scale_x_continuous(labels =label_number(suffix = deg_c)) +labs(y ="Exponent of ISD",title ="Within area*season Relationship",x ="5-Year Rolling Average Bottom Temperature")btemp_local
A spring spectra is not any more/less steep if spring is hotter or colder than usual. However, Fall is consistently steeper than spring, and fall is always hotter than spring.
Landings Relationship
Since landings are annual there won’t be any seasonal interaction, only an intercept adjustment for the survey season.
Code
# Clean up the plot:# Plot the predictions over datamod2_preds <-as.data.frame(ggpredict( mod2, terms =~ total_weight_lb * survey_area * season)) %>%mutate(survey_area =factor(group, levels = area_levels_long_all),season =factor(facet, levels =c("Spring", "Fall")))#### Trend Posthoc ####trend_df <-emtrends(object = mod2,~survey_area * season,var ="total_weight_lb")#summary(trend_df, infer= c(T,T))# Drop effect fits that are non-significant ###mod2_emtrend <-emtrends(object = mod2,specs =~ survey_area*season,var ="log(total_weight_lb)") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T),group =str_c(survey_area, "X", season))# Plot over observed datamod2_preds %>%left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%filter(non_zero) %>%ggplot() +geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha =0.1) +geom_line(aes(x, predicted, color = season), linewidth =1) +geom_point(data = wtb_model_df,aes(total_weight_lb, b, color = season),alpha =0.4,size =1) +facet_grid(survey_area~., scales ="free") +scale_color_gmri() +scale_x_continuous(transform ="log10", labels =label_log(base =10)) +labs(y ="Body Mass Spectra Slope (b)",title ="Wigley Species, Body Mass Spectra",x ="Total Landings (lb.)")
Landings + Temperature Modeling Considerations
With the old landings data we had trouble determining attribution from the models with both temperature and landings simultaneously due to strong collinearities. This may no longer be the case, lets look.
The bigger issue looks to be the relationship between survey season + region and bottom temperature. There are other phenological dynamics (spring plankton bloom, seasonal migrations, spawning, stratification etc.) that happen seasonally that I am hesitant to use temperature alone without controlling for season.
If we check variance inflation factors for the model we currently have nearly all predictors have very high VIF values, beyond what is the typical cutoff.
Which is unfortunate, because it seems like its a pretty reasonable model otherwise.
Code
# Drop effect fits that are non-significant #### summary(mod2)# summary(temp_landings_mod)check_model(mod2)
Everything below is digging into different nuances of the approach or the community composition. These things are of secondary importance and will be moved out.
Spiny Dogfish Migrations
Spiny dogfish constitute a majority fraction of biomass in the MAB during the Spring survey, but are subsequently absent during the fall season.
Their numbers have also been increasing over time and they are unquestionably the dominant “large” species caught in the survey.
Their growing abundance, particularly their spring abundance in the Mid-Atlantic Bight, acts to shallow the spectra as they are one of the larger species routinely sampled by the trawl survey.
If we remove them from the estimation and re-run the slopes this is what changes.
Once removed we see that the Mid-Atlantic Bight is no longer becoming less steep over time, with no long-term trend. Spectra in Georges Bank without the inclusion of spiny dogfish are now declining in both seasons.
Code
# # Run MAB with and without spiny dogfish# # Run the year, season, area fits for the filtered data# no_dogs_bmspectra <- group_binspecies_spectra(# ss_input = filter(wigley_in, comname != "spiny dogfish"),# grouping_vars = c("est_year", "season", "survey_area"),# abundance_vals = "numlen_adj",# length_vals = "length_cm",# use_weight = TRUE,# isd_xmin = 1,# global_min = TRUE,# isd_xmax = NULL,# global_max = FALSE,# bin_width = 1,# vdiff = 2)# # Save those# write_csv(# no_dogs_bmspectra,# here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv"))
Code
# Read them in, do some plottingno_dogs_bmspectra <-read_csv( here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv")) %>%mutate(survey_area =factor(survey_area, levels = area_levels))# Format a littleno_dogs_bmspectra <- no_dogs_bmspectra %>%mutate(est_year =as.numeric(est_year))# Join them togetherdfish_results <-bind_rows(list("Wigley Full"= wigley_bmspectra,"Wigley w/o Spiny Dogfish"= no_dogs_bmspectra), .id ="community") %>%mutate(metric ="bodymass_spectra") %>%mutate(survey_area =fct_relevel(survey_area, area_levels))# Get trends for no dogfishnodfish_mod <-lmer(formula = b ~ est_year * survey_area * season + (1|est_year),data = no_dogs_bmspectra)bmspectra_mod_wig <-lmer(formula = b ~ est_year * survey_area * season + (1|est_year),data = wigley_bmspectra)# Put predictionsdogfish_sens_predictions <-bind_rows(list("bodymass_spectra-Wigley Full"=get_preds_and_trendsignif(bmspectra_mod_wig),"bodymass_spectra-Wigley w/o Spiny Dogfish"=get_preds_and_trendsignif(nodfish_mod)), .id ="model_id") %>%separate(model_id, into =c("metric", "community"), sep ="-") %>%mutate(metric =fct_rev(metric))# Plot the differencesggplot() +geom_ribbon(data =filter(dogfish_sens_predictions, non_zero ==TRUE),aes(x, ymin = conf.low, ymax = conf.high, fill = season),linewidth =1, alpha =0.35) +geom_point(data = dfish_results,aes(est_year, b, color = season, alpha =I(if_else(survey_area %in%c("MAB", "GB"), 1, 0.2))),size =0.8) +geom_line(data =filter(dogfish_sens_predictions, non_zero ==TRUE),aes(x, predicted, color = season,alpha =I(if_else(survey_area %in%c("MAB", "GB"), 1, 0.2))),linewidth =1) +scale_fill_gmri() +scale_color_gmri() +facet_grid(survey_area ~ metric*community, scales ="free") +labs(x ="Year",y ="Exponent of Size Spectra",title ="Shift in Spectra Trend(s) with Exclusion of Spiny Dogfish",subtitle ="Spring = Dogfish in MAB, Fall = Dogfish in GB")
Spiny Dogfish Seasonal Dstributions
quicl GAM check to see if the distribution shifts seasonally.
Story Thoughts
Below is a table of some papers operating in the same area or around the same topic that are particularly relevant:
Code
shelf_papers <-tribble(~"Author", ~"Year", ~"Conjecture","Friedland et al.", 2024, "Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.","Krumsick", 2020, "Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios")shelf_papers %>%gt()
Author
Year
Conjecture
Friedland et al.
2024
Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.
Krumsick
2020
Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios
There are a couple core ideas about what changes have occurred or have been ocurring over the Northeast Shelf with relevance to the community that is sampled by the trawl survey: